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We present two new results regarding damage spreading in ferromagnetic Ising models. First, we 
show that a damage spreading transition can occur in an Ising chain that evolves in contact with a 
thermal reservoir. Damage heals at low temperature and spreads at high T. The dynamic rules for 
the system's evolution for which such a transition is observed are as legitimate as the conventional 
rules (Glauber, Metropolis, heat bath). Our second result is that such transitions are not always in 
the directed percolation universality class. 
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I. INTRODUCTION 



A system is said to exhibit damage spreading (DS) if 
the "distance" between two of its replicas, that evolve 
under the same thermal noise but from slightly different 
initial conditions, increases with time. Even though DS 
was first introduced in the context of biologically moti- 
vated dynamical systems Q, it has evolved into an im- 
portant tool in physics. It is used in equilibrium for 
measuring accurately dynamic exponents and also out 
of equilibrium, to study the influence of initial condi- 
tions on the temporal evolution of various systems. In 
particular, one hoped that DS could be used to identify 
"phases" of chaotic behavior in systems with no intrin- 
sic dynamics, such as Ising ferromagnets fj|,|| and spin- 
glasses Such hopes were dampened when it was re- 
alized that different algorithmic implementations of the 
same physical system's dynamics (such as Glauber ver- 
sus heat bath or Metropolis Monte Carlo) can have dif- 
ferent DS properties This implies that DS is not 
an intrinsic property of a system ||, since two equally 
legitimate algorithms yield contradictory results. This 
problem was addressed recently in Q , where we realized 
that one can define "phases" on the basis of their DS 
properties in an algorithm-independent manner. To do 
this one must, however, consider simultaneously the en- 
tire set A of possible algorithms (dynamic procedures) 
that are consistent with the physics of the model studied 
(such as detailed balance, interaction range and symme- 
tries). Every system must belong to one of three possible 
DS phases, depending on whether damage spreads for all, 
none or a part of the members of the set A. 

Once we have been led to consider a large family of al- 
gorithms, it was natural to revisit an old question, such as 
the possibility for DS in the one-dimensional (1-d) Ising 
ferromagnet. In this case all conventional dynamic pro- 
cedures agree that damage does not spread. We show 
here that once the family of dynamic procedures is ex- 
tended in the spirit explained above, a DS transition is 
possible in the 1-d Ising model. Having found such a 
DS transition, it is again natural to investigate to which 



universality class it belongs. So far this issue could be 
addressed only for the 2-d case; since it is much easier to 
obtain high-quality numerical data in 1-d, we were able to 
test carefully a conjecture of Grassberger ||, according 
to which the generic universality class of damage spread- 
ing transitions is directed percolation (DP). This indeed 
is correct, but we discovered that if the dynamics that 
is being used has certain symmetries, the DS transition 
is not in the DP class. Interestingly this is the case for 
Glauber dynamics of the H = Ising model, for which 
the DS transition is non-DP. 

We start by reviewing briefly j|,|7| the conventional al- 
gorithms - Glauber, heat bath (HB) and Metropolis - 
and show that they form a particular subset of some gen- 
eral set of legitimate rules A. All members of A satisfy 
detailed balance with respect to the same Hamiltonian; 
hence all these rules generate the same equilibrium en- 
semble as the conventional algorithms and are equally 
legitimate to mimic the temporal evolution of an Ising 
system in contact with a thermal reservoir. Next, we in- 
troduce two "new" dynamic rules, which constitute just 
another subset of A, and show that for these two rules a 
DS transition does occur in the l-d Ising model. More- 
over, as we show in the example of the second rule, an 
additional Z2 symmetry of the DS order parameter leads 
to a transition that is not in the DP universality class. 



II. PREVIOUS WORK, WITH CONVENTIONAL 
ALGORITHMS 



Denote the site which is being updated by i and the 
set of its neighbors by j. The energy at time t is given 
by 



H 



(1) 



where = J/ksT and o~i(t) = ±1. Define a transition 
probability Pi(t): 



Pi(t) 



ohi(t) 



e hi(t) _|_ e -hi(t) ■ 



(2) 



The update rules of HB, Glauber and Metropolis dynam- 
ics are expressed in terms of random numbers z = zi (t) , 
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selected with equal probability from the interval [0,1]. 
The rule for standard HB is 



<Ti(t + 1) = signet) - z]. 



(3) 



A different dynamic process is obtained by generating at 
each site two independent random numbers, z + and z_, 
and using the first if <Ti(t) = +1 and the second when 
crj(t) = — 1. The rules of this uncorrelated HB dynamics 
may be written as 



■H = / sign[pi(t) 
sign[pi(f) 



z_l_] if o"i(t) = +1 

Z_] if (7,(0 = -1 



(4) 



Glauber dynamics uses only one random number per site: 



<Ti(t+l) = 



-sign[pi(t) - z] 
-sign[l -Pi(t) 



if 

if cri(t) 



(5) 



This rule can be expressed in the form of (||) but with 
the two random numbers completely anticorrelated, i.e., 
z+ + Z— = 1. 

Finally the rules for Metropolis dynamics read 



Oi(t+l) 



-sign[p^(t) 
-sign[Pi (t) 



if o-i(i) 
if (Tilt) 



(6) 



where (i) = min(l, e T2 ' li(t) ). 

It is easy to show that given (T i _ 1 (i), 0"i(t), (Tj +1 (t), the 
probability to get <7j (t + 1) = +1 is the same for standard 
HB, uncorrelated HB, and Glauber dynamics^]. Hence, 
by observing the temporal evolution of a single Ising sys- 
tem, one cannot tell by which of these methods was its 
trajectory in configuration space generated. 
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TABLE I. Two-point correlations in the one-dimensional 
Ising model for various dynamic rules. We used the notation 
k = tanh ■ 
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The difference between these dynamics may become 
evident only when we observe the evolution of two repli- 
cas, i.e., study damage spreading! Indeed Stanley et al 
H and also Mariz et al || found, using Glauber dy- 
namics, that damage spreads for the 2-d Ising model for 
T > T c ; similarly for Metropolis dynamics Q. More re- 
cently Grassberger |[o| claimed that the DS transition 
occurs slightly below T c for Glauber which was also ob- 
served in the corresponding mean field theory On 
the other hand damage does not spread at any temper- 
ature with standard HB for neither the 2-d Q nor the 
3-d Ising models ||. The 3-d model did exhibit DS for 
T > T* with T* < T c when Metropolis |2) and Glauber 
dynamics were used. In the l-d Ising model with 
HB, Glauber or Metropolis dynamics no damage spread- 
ing has been observed. 



III. GENERAL CLASS OF DYNAMIC 
PROCEDURES FOR THE ISING MODEL 

The dynamic rules considered here for the l-d Ising 
model consist of local updates, where a random variable 
r = ±1 is assigned to the spin tr^: 



a t {t + 1) := r <Ti _ 1 ( t ) i<Ti ( t ) j(74+1 ( t ) 



(7) 



This random variable is generated in some probabilistic 
procedure using one or several random numbers. Like 
in the conventional algorithms discussed above, we allow 
the random variable to depend only on the values taken 
at time t by the updated spin itself and the spins with 
which it interacts (i.e., its nearest neighbors). The set 
of all one-point functions (r ai _ 1 ,a z ,a i+1 ) determines the 
transfer matrix of a single system. Here (...) denotes 
the average over many independent realizations of ran- 
dom numbers. The simultaneous evolution (and, hence, 
DS) of two replicas {a} and {a'} is, however, governed 
by a joint transfer matrix of the two systems which, in 
turn, is completely determined by the two-point func- 
tions (r ff4 _ 1)(7 - 4)<Ti+I JV ff ( )ff < ). In general, n-point func- 
tions determine the joint transfer matrix of n replicas. 
An important requirement is that all correlation func- 
tions have to be invariant under the symmetries of the 
model |J. For a homogeneous Ising chain in zero field 
these symmetries are invariance under reflection 



1 <Ti_i ,a i ,n i + 1 1 a 



V<Ti-l,tTi,tTi+l) — ( 7 V i + i,£r i ,<T i _i) , 



and global inversion of all spins [Zi symmetry): 

(?Vi_i,CTi,CTi + i) — ~ ( r -o- i _ 1 ,-cr i ,-cr i+1 ) ] 



ai,a i+1 1 u' 



r' 



(8) 



(9) 



*For Metropolis dynamics the transition probabilities are 
different 



For both HB and for Glauber dynamics the one-point 
functions are given by 
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_ l )=2p i -l. 



(10) 



The corresponding transfer matrices for single systems 
are, hence, identical. On the other hand the two-point 
functions for HB and Glauber dynamics are different so 
that damage evolves differently (see Table I). Still, dam- 
age does not spread in l-d for any of these algorithms at 
any temperature. 



IV. DYNAMIC RULE FOR WHICH DAMAGE 
DOES SPREAD IN 1-D 

Consider the following dynamics for the l-d Ising 
model: 



+sign(p 4 - z) if a t 
-sign(l - pi - z) if Oi 



1 — Ci+l 
-1 ^ O-'H-l 



(11) 



As can be checked easily, this dynamical rule yields the 
same one-point correlations as in eq. (|l0|). Therefore, 
the evolution of a single replica using this rule cannot 
be distinguished from that of Glauber or HB dynamics. 
However, the two-point correlations (and therewith dam- 
age spreading properties) are different (see Table I). 

Unlike Glauber and HB, this dynamics does exhibit a 
damage spreading transition in l-d. This can be seen as 
follows. At T — oo eq. (11) reduces to 



r<j,_ 1 ,<j u a l+1 = <7j_i<7j+i sign(i - z) . 



(12) 



which implies that the local damage Ai(t) — 1— 5 ai ^ t y a i^ 
evolves deterministically: 



Aj(t + 1) 



if A 4 _i(t) = Ai+i(t) 

1 if Aj_i (i) ^ A i+1 (i) 



(13) 



Since this is exactly the update rule of a Domany-Kinzel 
model (lj in the active phase (with pi = 1 and p 2 = 0), 
we conclude that for T = oo damage spreads. On the 
other hand, for T — eq. ( |TT| ) reduces to 



sign(z 



if (Jj_i = cr,;+i 

|) if C7j_l ^ <T i+ l 



(14) 



In this case damage evolves probabilistically and cannot 
be viewed as an independent process. One can, however, 
show that the expectation value to get damage at site i, 
averaged over many realizations of random numbers sat- 
isfies the inequality (A,-(t + 1)) < |(Aj_i(i) + A i+1 (t)), 
that is (A(t + 1)) < (A(i)). This means that for T = 
damage does not spread. In fact, simulating the spread- 
ing process one observes a DS transition at finite temper- 
ature. A typical temporal evolution near the transition 
is shown in Fig. ^. 

In order to determine the critical exponents that char- 
acterize the DS transition, we perform dynamic Monte- 
Carlo simulations jlq]. Two replicas are started from 





replica 1 replica 2 damage 

FIG. 1. Temporal evolution of damage in the l-d Ising 
model of size 200 with the dynamics of eq. (11) near the DS 
transition J/fcsT*=0.2305. Each configuration is represented 
by a row of pixels and time goes downwards. The two replicas 
are started from identical initial conditions. At an early time, 
a damage of 5 sites is inserted in the center. 
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FIG. 2. Numerical results for the l-d Ising model with 
A: the dynamics of eq. (11) and B: the dynamics of eq. (18). 
The measured quantities are explained in the text. 

identical random initial conditions, where one damaged 
site is inserted at the center. Both replicas then evolve 
according to the dynamic rules of the system using the 
same set of random numbers. In order to minimize finite- 
size effects, we simulate a large system of 5000 sites with 
periodic boundary conditions. For various temperatures 
we perform 10 6 independent runs up to 1500 time steps. 
However, in many runs damage heals very soon so that 
the run can be stopped earlier. As usual in this type 
of simulations, we measure the survival probability P(t), 
the number of damaged sites A(i), and the mean-square 
spreading of damage from the center R 2 (t) averaged over 
the active runs. At the DS transition, these quantities are 
expected to scale algebraically in the large time limit: 

-5 



p(t) ~ r 



A(t) - f> , 



R 2 (t)~t z . (15) 



The critical exponents 5,r),z are related to the den- 
sity exponent /3 and the scaling exponents v\\ by 
5 = f3/v\\, z — 2i^/V|| and obey the hyperscaling re- 
lation AS + 2i] = dz. At criticality, the quantities ( |T5| ) 
show straight lines in double logarithmic plots. Off crit- 
icality, these lines are curved. Using this criterion we 
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estimate the critical temperature for the DS transition 
by J/k B T* = 0.2305(5). The exponents 5, 77, and z 
are measured at criticality while the density exponent 
[3 is determined off criticality by measuring the station- 
ary Hamming distance A(T) ~ (T — T* Y 3 in the spread- 
ing phase. The results of our simulations are shown in 
Fig. g. From the slopes in the double logarithmic plots 
we obtain the estimates 8 — 0.165(5), r\ = 0.315(10), 
z = 1.29(3), and (3 = 0.26(2) which are in fair agreement 
with the known []l6| exponents for directed percolation 
5 = 0.15947(3), T] = 0.31368(4), z = 1.26523(4), and 
j3 = 0.27649(4). We therefore conclude that in agree- 
ment with Grassberger's conjecture [gj, the DS transition 
belongs to the DP universality class. This is very plau- 
sible; as far as the damage variable is concerned there 
is a single absorbing state (of no damage at all) and the 
transition is from a phase in which the system ends up in 
this state to one in which it does not, just as is the case 
for DP. 



V. DAMAGE SPREADING TRANSITION WITH 
NON-DP EXPONENTS 



Different critical properties are expected Jl7|-22[ for 
rules with two distinct absorbing states (of the damage 
variables!) related by symmetry. It is important to note 
that the Z 2 symmetry of the Ising system does not suffice 
- inverting all spins in both replicas does not change the 
damage variable (the Hamming distance between the two 
configurations). Therefore, we are looking for dynamic 
rules which (a) have two types of absorbing states; one 
with no damage and the other with full damage. Fur- 
thermore, (b) the two play completely symmetric roles. 
One can see that both (a) and (b) hold for rules that 
satisfy the condition 



(16) 



The immediate consequence of this condition is that 
if a configuration {o~(t)} evolves in one time step into 
{er(< + 1)}, then the spin-reversed configuration {— o~(t)} 
will evolve into precisely {— a(t+l)}. Imagine now simul- 
taneous evolution of two replicas with initial states {a} 
and {a 1 }, giving rise to a damage field {A}. Reversal 
of the initial state on one of the replicas will give sign- 
reversed spin states on this replica and hence the damage 
field {—A} will evolve. Thus, for rules that satisfy con- 
dition (16), the damage variable has an Zi symmetry. A 
particular consequence of this symmetry is that if two ini- 
tial states are the exact sign-reversed of one another, this 
will persist at all subsequent times. Therefore, inasmuch 
as A = (no damage) is an absorbing state, so is the 
situation of full damage, A = 1. For systems with such 
symmetry we expect the DS transition (if it exists) to 
exhibit non-DP behavior. 





(a) (b) (c) 

FIG. 3. ^-symmetric damage spreading transition. Two 
replicas with 200 sites are started from identical random ini- 
tial conditions. At an early time 5 damaged sites are intro- 
duced in the center. For fixed temperature J/ksT = 0.25 a 
typical temporal evolution of damage is shown for (a) Glauber 
dynamics A = 1, (b) near the transition A* = 0.82 and (c) in 
the spreading regime A = 0. Because of the symmetry, islands 
of damaged sites can heal only at the edges. 

It is quite remarkable to note that Glauber dynamics 
satisfies eq. (|l^)! The Z 2 -symmetry of damage in the 
1-d Glauber model is illustrated in Fig. ||a. One can see 
that compact islands of damaged sites are formed because 
damage does not heal spontaneously inside such islands 
but only at the edges. However, as mentioned earlier, 
there is no DS transition in the 1-d Glauber model. 

Consider now a different dynamic rule: 



-sign(p; - z) if fTi-ifTiO-j+i = 1 
-sign(l -pi - z) if er i _iCT i (7. i+ i = -1 

(17) 



For this rule, which also satisfies eq. (16), we observe in 
simulations that damage always spreads (see Fig. Bp). In 
order to generate an ^-symmetric DS transition in 1-d, 
we use a rule that interpolates between this and Glauber. 
This can be done by introducing a second parameter 
< A < 1 and 'switching' between Glauber dynamics 
and rule ([l7|) as follows: in each update an additional 
random number z is generated. If z > A, rule (17]) is ap- 
plied, otherwise Glauber dynamics is used. This mixed 
dynamics can be expressed as 



-sign (p, 
-sign(l 



-z) 

Pi - 



if y 
if y 



(18) 



where y = ia i [(l + cr i _ 1 CT i+1 ) + (l-cr l _ 1 (T l+1 )sign(A-z)]. 
Again this rule leads to the one-point correlations of eq. 
(|10|), i.e., the temporal evolution of a single replica is 
the same as in Glauber and HB dynamics. However, 
varying A (at fixed T) we find a critical value A* where 
a DS transition occurs. A typical temporal evolution of 
damage near the transition is shown in Fig. |3jo. 

Since 'damage' and 'no damage' play a symmetric role, 
the Hamming distance A (the density of damaged sites) 
cannot be used as an order parameter. Instead one has 
to use the density of kinks N (domain walls) between 
damaged and healed domains. By definition, the num- 
ber of kinks is conserved modulo two which establishes 
a parity conservation law. As can be seen in Fig. |^, 
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two processes compete with each other: kinks annihilate 
mutually (2X — > 0) and already existing kinks branch 
into an odd number of kinks (X — > 3X, 5X, . . .). Both 
processes resemble a branching annihilating walk with an 
even number of offspring. This branching process has a 
continuous phase transition that belongs to the so-called 
parity-conserving (PC) universality class. Phase tran- 
sitions of this type have been observed in a variety of 
models, including certain probabilistic cellular automata 
|l7f , nonequilibrium kinetic Ising models with combined 
zero- and infinite-temperature dynamics |18|, interact- 
ing monomer-dimer models fL9| , branching-annihilating 
random walks |2(|| a nd certain lattice models with two 
absorbing states |2l|. In all these models the symmetry 
appears either as a parity conservation law or as an ex- 
plicit ^-symmetry among different absorbing phases. A 
field theory describing PC transitions is currently devel- 
oped in K2j. 

The PC universality class is characterized by the ex- 
ponents S = 0.285(5), 77 = 0.00(1), z = 1.15(1), and 
(3 = 0.92(2). In fact, repeating the numerical simulations 
described above for J/k B T = 0.25 and A* = 0.82(1) 
(see Fig. ||), we obtain the estimates S = 0.295(10), 
7] = 0.01(2), z = 1.17(3), and j3 = 0.86(5), which are in 
fair agreement with the known values. We therefore con- 
clude that the DS transition observed for the dynamics 
of eq. ( fl8| ) belongs to the PC universality class. Fur- 
thermore, our findings imply that the DS transitions ob- 
served |l^] for the 2-d Ising model with Glauber dynamics 
should also exhibit PC exponents (remember: d = 2) in 
zero field, and cross over to (2-d) DP values when a field 
is switched on. 



[12] U. M. S. Costa, J. Phys. A 20, L 583 (1987). 
[13] G. La Caer, J. Phys. A 22, L 647 (1989); Physica A 159, 
329 (1989). 

[14] E. Domany and W. Kinzel, Phys. Rev. Lett. 53, 447 
(1984). 

[15] P. Grassberger and A. de la Torre, Ann. Phys. (N.Y.) 

122, 373 (1979). 
[16] I. Jensen, Phys. Rev. Lett. 77, 4988 (1996). 
[17] P. Grassberger, F. Krause, and T. von der Twer, J. Phys. 

A 17 L105 (1984); P. Grassberger, J. Phys. A 22, L1103 

(1984). 

[18] N. Menyhard, J. Phys. A 27, 6139 (1994); B. Menyhard 
and G. Odor, J. Phys. A 29, 7739 (1996). The kinetic 
Ising models introduced in these papers have mixed zero- 
and infinite-temperature dynamics and are different from 
the usual Ising model at finite temperature discussed in 
the present work. 

[19] H. H. Kim and H. Park, Phys. Rev. Lett. 73, 2579 (1994); 

H. Park, M. H. Kim, and H. Park, Phys. Rev. E 52, 5664 

(1995) . 

[20] I. Jensen, J. Phys. A 26 3921 (1993); D. ben-Avraham, 
F. Leyvraz, and S. Redner, Phys. Rev. E 50, 1843 (1994); 

I. Jensen, Phys. Rev. E 50 3623 (1994); D. Zhong and 
D. ben-Avraham, Phys. Lett. A 209, 333 (1995). 

[21] H. Hinrichsen, Phys. Rev. E 55, 219 (1997). 

[22] J. Cardy and U. C. Tauber, Phys. Rev. Lett. 77, 4780 

(1996) . 



We thank D. Stauffer for sharing with us his knowledge 
of the DS literature and for encouragement. This work 
was supported by The Minerva Foundation and by the 
Germany-Israel Science Foundation (GIF). 



[1] S.A. Kauffman, J. Theor. Biol. 22, 437 (1969). 

[2] P. Grassberger, Physica A 214, 547 (1995). 

[3] M. Creutz, Ann. Phys. 167, 62 (1986). 

[4] H. Stanley, D. Stauffer, J. Kertesz and H. Herrmann, 

Phys. Rev. Lett. 59, 2326 (1987). 
[5] B. Derrida and G. Weisbuch, Europhys. Lett. 4, 657 

(1987). 

[6] A. M. Mariz, H. J. Herrmann and L. de Arcangelis, J. 

Stat. Phys. 59, 1043 (1990). 
[7] N. Jan and L. de Arcangelis, Ann. Rev. Comp. Phys. 1,1 

(ed. D. Stauffer, World Scientific, Singapore 1994). 
[8] P. Grassberger, J. Stat. Phys. 79, 13 (1995). 



[9] H. Hinrichsen , S. Weitz, and E. Domany, preprint cond- 
mat/9611085| , submitted to J. Stat. Phys. 
[10] P. Grassberger, J. Phys. A 28, L 67 ( 1995). 
[11] T. Vojta, preprint |co~nd-mat/9610084| (1996). 



5 



